Clustering Genes

K-means Clustering

set.seed(123)
wss = sapply(1:15, function(k) kmeans(datExpr_redDim, k, iter.max=100, nstart=25,
                                      algorithm='MacQueen')$tot.withinss)
plot(wss, type='b', main='K-Means Clustering')
best_k = 4
abline(v = best_k, col='blue')

datExpr_k_means = kmeans(datExpr_redDim, best_k, iter.max=100, nstart=25)
clusterings[['km']] = datExpr_k_means$cluster

Hierarchical Clustering

Chose k=6 as best number of clusters.

Clusters seem to be able to separate ASD and control samples pretty well and there are no noticeable patterns regarding sex, age or brain region in any cluster.

Younger ASD samples seem to be more similar to control samples than older ASD samples (pink cluster has most of the youngest samples). The yellow cluster is made of young ASD samples.

h_clusts = datExpr_redDim %>% dist %>% hclust
plot(h_clusts, hang = -1, cex = 0.6, labels=FALSE)

best_k = 4
clusterings[['hc']] = cutree(h_clusts, best_k)

Consensus Clustering

Samples are grouped into two big clusters, one small, one of only six genes and two outliers, and two outliers, the first big cluster has three main subclusters, two small subclusters and two outliers, and the second one two main subclustersand five small subclusters.

*Output plots in clustering_genes_03_20 folder

Independent Component Analysis

Following this paper’s guidelines:

0. Run PCA and keep enough components to explain 60% of the variance

1. Run ICA with that same number of nbComp as principal components kept to then filter them

2. Select components with kurtosis > 3

3. Assign genes to clusters with FDR<0.01 using the fdrtool package

ICA_output = datExpr_redDim %>% runICA(nbComp=ncol(datExpr_redDim), method='JADE')
signals_w_kurtosis = ICA_output$S %>% data.frame %>% dplyr::select(names(apply(ICA_output$S, 2, kurtosis)>3))
ICA_clusters = apply(signals_w_kurtosis, 2, function(x) fdrtool(x, plot=F)$qval<0.01) %>% data.frame
ICA_clusters = ICA_clusters[colSums(ICA_clusters)>1]

clusterings[['ICA_min']] = rep(NA, nrow(ICA_clusters))
ICA_clusters_names = colnames(ICA_clusters)
for(c in ICA_clusters_names) clusterings[['ICA_min']][ICA_clusters[,c]] = c

Leaves most of the observations (~80%) without a cluster (same result as when considering the first PC):

ICA_clusters %>% rowSums %>% table

## .
##    0    1    2    3    4    5 
## 2350  420  117   32    8    1

ICA_clusters %>% mutate(cl_sum=rowSums(.)) %>% as.matrix %>% melt %>% ggplot(aes(Var2,Var1)) + 
  geom_tile(aes(fill=value)) + xlab('Clusters') + ylab('Samples') + 
  theme_minimal() + theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + coord_flip()

WGCNA

best_power=21

best_power = datExpr_redDim %>% t %>% pickSoftThreshold(powerVector = seq(1, 30, by=1))

##    Power SFT.R.sq   slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.7890  1.7600          0.818 1360.00   1450.00 1880.0
## 2      2   0.3970  0.4290          0.486  817.00    871.00 1390.0
## 3      3   0.0295 -0.0815          0.133  549.00    571.00 1090.0
## 4      4   0.3490 -0.3280          0.520  394.00    397.00  885.0
## 5      5   0.5340 -0.5000          0.660  295.00    285.00  735.0
## 6      6   0.6380 -0.6140          0.782  228.00    211.00  622.0
## 7      7   0.6630 -0.7190          0.810  181.00    161.00  534.0
## 8      8   0.7110 -0.7890          0.849  146.00    123.00  462.0
## 9      9   0.7450 -0.8490          0.893  120.00     97.00  404.0
## 10    10   0.7620 -0.8940          0.915   99.60     77.40  356.0
## 11    11   0.7750 -0.9340          0.927   83.80     62.40  316.0
## 12    12   0.7880 -0.9730          0.941   71.20     51.30  282.0
## 13    13   0.8050 -1.0000          0.950   61.00     42.40  253.0
## 14    14   0.8160 -1.0300          0.960   52.70     35.10  228.0
## 15    15   0.8240 -1.0600          0.964   45.90     29.80  207.0
## 16    16   0.8260 -1.1000          0.965   40.10     25.30  188.0
## 17    17   0.8300 -1.1300          0.966   35.30     21.70  172.0
## 18    18   0.8420 -1.1500          0.973   31.30     18.50  157.0
## 19    19   0.8390 -1.1800          0.970   27.80     16.00  144.0
## 20    20   0.8480 -1.1900          0.974   24.80     13.90  133.0
## 21    21   0.8580 -1.2100          0.979   22.20     12.20  123.0
## 22    22   0.8600 -1.2300          0.979   20.00     10.60  114.0
## 23    23   0.8650 -1.2400          0.981   18.00      9.35  105.0
## 24    24   0.8690 -1.2500          0.983   16.30      8.30   97.8
## 25    25   0.8730 -1.2600          0.987   14.80      7.38   91.1
## 26    26   0.8800 -1.2700          0.989   13.50      6.61   84.9
## 27    27   0.8810 -1.2800          0.990   12.30      5.91   79.4
## 28    28   0.8840 -1.2900          0.990   11.30      5.29   74.3
## 29    29   0.8840 -1.3000          0.990   10.40      4.75   69.6
## 30    30   0.8880 -1.3100          0.991    9.52      4.28   65.3

network = datExpr_redDim %>% t %>% blockwiseModules(power=best_power$powerEstimate, numericLabels=TRUE)
clusterings[['WGCNA']] = network$colors
names(clusterings[['WGCNA']]) = rownames(datExpr_redDim)

It leaves 484 genes without a cluster (even more than with the first PC)

clusterings[['WGCNA']] %>% table

## .
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17 
## 484 423 384 337 329 179 147 132 120  64  55  54  46  36  32  32  30  24 
##  18 
##  20

Gaussian Mixture Models with hard thresholding

Number of clusters that resemble more Gaussian mixtures = 32 but 25 clusters gets a very similar BIC value

n_clust = datExpr_redDim %>% Optimal_Clusters_GMM(max_clusters=50, criterion='BIC', plot_data=FALSE)
plot(n_clust, type='l', main='Bayesian Information Criterion to choose number of clusters')

best_k = 25
gmm = datExpr_redDim %>% GMM(best_k)
clusterings[['GMM']] = gmm$Log_likelihood %>% apply(1, which.max)

Plot of clusters with their centroids in gray

gmm_points = rbind(datExpr_redDim, setNames(data.frame(gmm$centroids), names(datExpr_redDim)))
gmm_labels = c(clusterings[['GMM']], rep(NA, best_k)) %>% as.factor
ggplotly(gmm_points %>% ggplot(aes(x=PC1, y=PC2, color=gmm_labels)) + geom_point() + theme_minimal())

Trying with 9 clusters

best_k = 9
gmm = datExpr_redDim %>% GMM(best_k)
clusterings[['GMM_2']] = gmm$Log_likelihood %>% apply(1, which.max)
clusterings[['GMM_2']] %>% table

## .
##   1   2   3   4   5   6   7   8   9 
## 305 239 220 345 295 281 350 475 418

Trying with 2 clusters to see if it can separate the two clouds of points (it can’t)

best_k = 2
gmm = datExpr_redDim %>% GMM(best_k)
clusterings[['GMM_3']] = gmm$Log_likelihood %>% apply(1, which.max)

gmm_points = rbind(datExpr_redDim, setNames(data.frame(gmm$centroids), names(datExpr_redDim)))
gmm_labels = c(clusterings[['GMM_3']], rep(NA, best_k)) %>% as.factor
ggplotly(gmm_points %>% ggplot(aes(x=PC1, y=PC2, color=gmm_labels)) + geom_point() + theme_minimal())

Manual clustering

Separate the two clouds of points by a straight line. There seems to be a difference in the mean expression of the genes between clusters but not in their standard deviation.

manual_clusters = as.factor(as.numeric(1.5*datExpr_redDim$PC1 + 0.5 > datExpr_redDim$PC2))
datExpr_redDim %>% ggplot(aes(x=PC1, y=PC2, color=manual_clusters)) + geom_point() + 
  geom_abline(slope=1.5, intercept=0.5, color='gray') + theme_minimal()

names(manual_clusters) = rownames(datExpr_redDim)

clusterings[['Manual']] = manual_clusters

Cluster 2 has a slightly bigger mean expression than cluster 1 and a smaller standard deviation.

manual_clusters_data = cbind(apply(datExpr_redDim, 1, mean), apply(datExpr_redDim, 1, sd), 
                             manual_clusters) %>% data.frame
colnames(manual_clusters_data) = c('mean','sd','cluster')
manual_clusters_data = manual_clusters_data %>% mutate('cluster'=as.factor(cluster))
p1 = manual_clusters_data %>% ggplot(aes(x=mean, color=cluster, fill=cluster)) + 
  geom_density(alpha=0.4) + theme_minimal()
p2 = manual_clusters_data %>% ggplot(aes(x=sd, color=cluster, fill=cluster)) + 
  geom_density(alpha=0.4) + theme_minimal()
grid.arrange(p1, p2, ncol=2)

rm(wss, datExpr_k_means, h_clusts, cc_output, cc_output_c1, cc_output_c2, best_k, ICA_output, 
   ICA_clusters_names, signals_w_kurtosis, n_clust, gmm, gmm_points, gmm_labels, network, 
   best_power, c, manual_clusters, manual_clusters_data, c2_sd, gmm_c2_sd, p1, p2, 
   pca_data_projection)

2D plots of clusterings

• Simple methods seem to only partition the space in buckets using information from the first component

• All clusterings except K-means manage to separate the two clouds at least partially, but no one does a good job

• WGCNA creates clusters inverted between clouds

• Manual classification + GMM by standard deviation clusters make sense

create_2D_plot('km_clust')

create_2D_plot('hc_clust')

create_2D_plot('cc_l1_clust')

create_2D_plot('cc_clust')

create_2D_plot('ica_clust')

create_2D_plot('gmm_clust')

create_2D_plot('gmm_clust2')

create_2D_plot('gmm_clust3')

create_2D_plot('wgcna_clust')

3D plots

create_3D_plot('ica_clust')

create_3D_plot('gmm_clust')

create_3D_plot('gmm_clust3')

create_3D_plot('wgcna_clust')